ENSO analysis new

Packages

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ───────────────────────────────────────────────────────────── tidyverse 1.3.2 ──✔ tibble  3.2.1     ✔ stringr 1.5.0
✔ readr   2.1.3     ✔ forcats 1.0.0
✔ purrr   1.0.1     ── Conflicts ──────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ tidyr::extract() masks raster::extract()
✖ dplyr::filter()  masks stats::filter()
✖ dplyr::lag()     masks stats::lag()
✖ purrr::map()     masks maps::map()
✖ dplyr::select()  masks raster::select()

library(sf)

Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE

library(here)

here() starts at /Users/ramkripa

library(weed)
library(rnaturalearth)
library(rnaturalearthdata)

Attaching package: ‘rnaturalearthdata’

The following object is masked from ‘package:rnaturalearth’:

    countries110

Loading In the Data

emdat_old <- read_csv('~/Desktop/projects/emdat_proj/full_emdat_geocoded_finally.csv')

Warning: One or more parsing issues, call `problems()` on your data frame for details, e.g.:
  dat <- vroom(...)
  problems(dat)Rows: 74035 Columns: 51── Column specification ────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (27): Dis No, Seq, Disaster Group, Disaster Subgroup, Disaster Type, Disaster Subtype, Disaster Sub...
dbl  (18): Year, Dis Mag Value, Start Year, Start Month, Start Day, End Year, End Month, End Day, Total ...
lgl   (5): Glide, Aid Contribution, Reconstruction Costs ('000 US$), Admin1 Code, uncertain_location_spe...
time  (1): Local Time
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

emdat_new <- read_csv('~/Desktop/projects/emdat_proj/data/geocoded_21_23_emdat.csv')

Rows: 4987 Columns: 54── Column specification ────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (27): Dis No, Seq, Glide, Disaster Group, Disaster Subgroup, Disaster Type, Disaster Subtype, Disast...
dbl (22): Year, Dis Mag Value, Start Year, Start Month, Start Day, End Year, End Month, End Day, Total D...
lgl  (5): AID Contribution ('000 US$), Latitude, Longitude, Local Time, uncertain_location_specificity
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

How many disasters were located?

How many locations were located?

All types of disasters

Restrict to dry and wet hydrological disasters

dry_list <- c("Drought", "Extreme temperature", "Wildfire")
wet_list <- c("Flood", "Landslide", "Storm")

emdat_hydrological <- emdat %>%
  filter((`Disaster Type` %in% dry_list) || (`Disaster Type` %in% wet_list)) %>%
  mutate(dry = (`Disaster Type` %in% dry_list), wet = (`Disaster Type` %in% wet_list))

Warning: There were 2 warnings in `filter()`.
The first warning was:
ℹ In argument: `(`Disaster Type` %in% dry_list) || ...`.
Caused by warning in `` (`Disaster Type` %in% dry_list) || (`Disaster Type` %in% wet_list) ``:
! 'length(x) = 74035 > 1' in coercion to 'logical(1)'
ℹ Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 1 remaining warning.

ENSO data

enso_data <- read_csv('~/Desktop/projects/emdat_proj/data/enso_data_copy2.csv')

Rows: 812 Columns: 7── Column specification ────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (4): MonthTxt, PeriodTxt, PeriodNum, Event
dbl (3): Year, MonthNum, Value
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

enso_data %>%
  group_by(Event) %>%
  count()

emdat_hydrological_enso %>%
  filter(!is.na(enso_Value)) %>%
  filter(!is.na(wet)) %>%
  mutate(enso_positive = (enso_Value > 0)) %>%
  group_by(`Dis No`) %>%
  summarize(enso_pos = any(enso_positive),enso_val = mean(enso_Value), is_wet = any(wet)) %>%
  ggplot(mapping = aes(x = enso_val, fill = is_wet)) +
    geom_histogram(stat = "density") +
    facet_grid(~is_wet) +
  ggtitle("Counts of Disasters compared to value of ONI with wet disasters in Blue and Dry in Red")

Warning: Ignoring unknown parameters: `binwidth`, `bins`, and `pad`

resy <- 2.5

emdat_hydrological_enso %>%
  left_join(event_num_months, by = c("Event")) %>%
  filter(!is.na(enso_Value)) %>%
  filter(!is.na(wet)) %>%
  filter(!is.na(lat)) %>%
  mutate(wet_factor = as_factor(wet)) %>%
  mutate(Event_factor = fct_recode(as_factor(Event), no_event = "N", negative = "-", positive = "+")) %>%
  mutate(wet_string = fct_recode(wet_factor, wet_disasters = "TRUE", dry_disasters = "FALSE")) %>%
  mutate(lat_box = (as.integer(lat) %/% resy) * resy , lng_box = (as.integer(lng) %/% resy) * resy) %>%
  mutate(dis_no_factor = as_factor(`Dis No`)) %>%
  group_by(lat_box, lng_box, Event_factor, wet_string, dis_no_factor) %>%
  summarize(num_months = first(num_months), num_disaster_month = mean(num_months_disaster)) %>%
  group_by(lat_box, lng_box, Event_factor, wet_string) %>%
  summarize(num_months = first(num_months), num_disaster_months = sum(num_disaster_month)) %>%
  mutate(disaster_rate = num_disaster_months / num_months) %>%
  pivot_wider(id_cols = c(lat_box, lng_box), names_from = c(Event_factor, wet_string), values_from = c(disaster_rate), values_fill = 0) %>%
  mutate(positive_wet_relative_risk = positive_wet_disasters / no_event_wet_disasters, 
         negative_wet_relative_risk = negative_wet_disasters / no_event_wet_disasters, 
         positive_dry_relative_risk =  positive_dry_disasters / no_event_dry_disasters,
         negative_dry_relative_risk = negative_dry_disasters / no_event_dry_disasters) %>%
  select(lat_box, lng_box, positive_wet_relative_risk, negative_wet_relative_risk, positive_dry_relative_risk, negative_dry_relative_risk) %>%
  pivot_longer(cols = c(positive_wet_relative_risk, negative_wet_relative_risk, positive_dry_relative_risk, negative_dry_relative_risk), names_to = 'risk', values_to = 'value') %>%
  mutate(risk_factor = as_factor(risk), values_cleaned = if_else(is.infinite(value), 18, value)) %>%
  mutate(relative_risk_cleaned = if_else(is.nan(values_cleaned), 1, values_cleaned)) %>%
  mutate(relative_risk_cleaned2 = if_else((relative_risk_cleaned) > 2, 2, relative_risk_cleaned)) %>%
  ggplot() +
  geom_sf(data = world_map) +
  geom_tile(mapping = aes(x = lng_box, y = lat_box, fill = relative_risk_cleaned2, width = resy)) +
  scale_fill_gradient2(midpoint = 1, limits = c(0, 2), low = 'cyan', high = 'red') +
  facet_grid(rows = vars(risk_factor))+
  ggtitle("Relative risk plot using calculations of the form [{num_disasters_positive_wet / num_months_positive} / {num_disasters_no_event_wet / num_months_no_event}]")

`summarise()` has grouped output by 'lat_box', 'lng_box', 'Event_factor', 'wet_string'. You can override using the `.groups` argument.`summarise()` has grouped output by 'lat_box', 'lng_box', 'Event_factor'. You can override using the `.groups` argument.

resy <- 2.5

emdat_hydrological_enso %>%
  left_join(event_num_months, by = c("Event")) %>%
  filter(!is.na(enso_Value)) %>%
  filter(!is.na(wet)) %>%
  filter(!is.na(lat)) %>%
  mutate(wet_factor = as_factor(wet)) %>%
  mutate(Event_factor = fct_recode(as_factor(Event), no_event = "N", negative = "-", positive = "+")) %>%
  mutate(wet_string = fct_recode(wet_factor, wet_disasters = "TRUE", dry_disasters = "FALSE")) %>%
  mutate(lat_box = (as.integer(lat) %/% resy) * resy , lng_box = (as.integer(lng) %/% resy) * resy) %>%
  mutate(dis_no_factor = as_factor(`Dis No`)) %>%
  group_by(lat_box, lng_box, Event_factor, wet_string, dis_no_factor) %>%
  summarize(num_months = first(num_months), num_disaster_month = mean(num_months_disaster)) %>%
  group_by(lat_box, lng_box, Event_factor, wet_string) %>%
  summarize(num_months = first(num_months), num_disaster_months = sum(num_disaster_month)) %>%
  mutate(disaster_rate = num_disaster_months / num_months) %>%
  pivot_wider(id_cols = c(lat_box, lng_box), names_from = c(Event_factor, wet_string), values_from = c(disaster_rate), values_fill = 0) %>%
  mutate(positive_wet_percent = positive_wet_disasters / (no_event_wet_disasters + positive_wet_disasters + negative_wet_disasters) , 
         negative_wet_percent = negative_wet_disasters / (no_event_wet_disasters + positive_wet_disasters + negative_wet_disasters), 
         positive_dry_percent =  positive_dry_disasters / (no_event_dry_disasters + positive_dry_disasters + negative_dry_disasters),
         negative_dry_percent =  negative_dry_disasters / (no_event_dry_disasters + positive_dry_disasters + negative_dry_disasters)
         ) %>%
  select(lat_box, lng_box, positive_wet_percent, negative_wet_percent, positive_dry_percent, negative_dry_percent) %>%
  pivot_longer(cols = c(positive_wet_percent, negative_wet_percent, positive_dry_percent, negative_dry_percent), names_to = 'typey', values_to = 'value') %>%
  mutate(type_factor = as_factor(typey), values = if_else(is.nan(value), 0, value)) %>%
  ggplot() +
  geom_sf(data = world_map) +
  geom_tile(mapping = aes(x = lng_box, y = lat_box, fill = values, width = resy)) +
  scale_fill_gradient(low = 'white', high = 'red') +
  facet_grid(rows = vars(type_factor))+
  ggtitle("Disaster Percent plot using calculations of the form [positive_wet / all_wet]")

`summarise()` has grouped output by 'lat_box', 'lng_box', 'Event_factor', 'wet_string'. You can override using the `.groups` argument.`summarise()` has grouped output by 'lat_box', 'lng_box', 'Event_factor'. You can override using the `.groups` argument.

# SPLIT YEAR INTO 2
# Apr to Sep
# Oct to Mar
## calculate relative risk within each box
### create boxes acc to resy X resy

resy <- 2.5

event_half_num_months <- enso_data %>%
  mutate(spring_summery = as_factor((MonthNum >= 4) & (MonthNum <= 9))) %>%
  mutate(half_year = fct_recode(spring_summery, spring_summer = "TRUE", fall_winter = "FALSE")) %>%
  group_by(Event, half_year) %>%
  summarize(num_months = n())

`summarise()` has grouped output by 'Event'. You can override using the `.groups` argument.

emdat_hydrological_enso %>%
  mutate(spring_summery = as_factor((end_month >= 4) & (end_month <= 9))) %>%
  mutate(half_year = fct_recode(spring_summery, spring_summer = "TRUE", fall_winter = "FALSE")) %>%
  left_join(event_half_num_months, by = c("Event", "half_year")) %>%
  filter(!is.na(enso_Value)) %>%
  filter(!is.na(wet)) %>%
  filter(!is.na(lat)) %>%
  mutate(wet_factor = as_factor(wet)) %>%
  mutate(Event_factor = fct_recode(as_factor(Event), no_event = "N", negative = "-", positive = "+")) %>%
  mutate(wet_string = fct_recode(wet_factor, wet_disasters = "TRUE", dry_disasters = "FALSE")) %>%
  mutate(lat_box = (as.integer(lat) %/% resy) * resy , lng_box = (as.integer(lng) %/% resy) * resy) %>%
  mutate(dis_no_factor = as_factor(`Dis No`)) %>%
  group_by(lat_box, lng_box, Event_factor, wet_string, dis_no_factor, half_year) %>%
  summarize(num_months = first(num_months), num_disaster_month = mean(num_months_disaster)) %>%
  group_by(lat_box, lng_box, Event_factor, wet_string, half_year) %>%
  summarize(num_months = first(num_months), num_disaster_months = sum(num_disaster_month)) %>%
  mutate(disaster_rate = num_disaster_months / num_months) %>%
  pivot_wider(id_cols = c(lat_box, lng_box, half_year), names_from = c(Event_factor, wet_string), values_from = c(disaster_rate), values_fill = 0) %>%
  mutate(positive_wet_relative_risk = positive_wet_disasters / no_event_wet_disasters, 
         negative_wet_relative_risk = negative_wet_disasters / no_event_wet_disasters, 
         positive_dry_relative_risk =  positive_dry_disasters / no_event_dry_disasters,
         negative_dry_relative_risk = negative_dry_disasters / no_event_dry_disasters) %>%
  select(lat_box, lng_box, half_year, positive_wet_relative_risk, negative_wet_relative_risk, positive_dry_relative_risk, negative_dry_relative_risk) %>%
  pivot_longer(cols = c(positive_wet_relative_risk, negative_wet_relative_risk, positive_dry_relative_risk, negative_dry_relative_risk), names_to = 'risk', values_to = 'value') %>%
  mutate(risk_factor = as_factor(risk), values_cleaned = if_else(is.infinite(value), 18, value)) %>%
  mutate(relative_risk_cleaned = if_else(is.nan(values_cleaned), 1, values_cleaned)) %>%
  mutate(relative_risk_cleaned_truncatedat2 = if_else((relative_risk_cleaned) > 2, 2, relative_risk_cleaned)) %>%
  select(lat_box, lng_box, half_year, risk_factor, relative_risk_cleaned, relative_risk_cleaned_truncatedat2) %>%
  ggplot() +
  geom_sf(data = world_map) +
  geom_tile(mapping = aes(x = lng_box, y = lat_box, fill = relative_risk_cleaned_truncatedat2, width = resy)) +
  scale_fill_gradient2(midpoint = 1, limits = c(0, 2), low = 'cyan', high = 'red') +
  facet_grid(rows = vars(risk_factor), cols = vars(half_year))+
  ggtitle("Half year Relative risk plot using calculations of the form [{num_disasters_positive_wet / num_months_positive} / {num_disasters_no_event_wet / num_months_no_event}]")

`summarise()` has grouped output by 'lat_box', 'lng_box', 'Event_factor', 'wet_string', 'dis_no_factor'. You can override using the `.groups` argument.`summarise()` has grouped output by 'lat_box', 'lng_box', 'Event_factor', 'wet_string'. You can override using the `.groups` argument.

NA

Dry Disasters

resy <- 2.5

emdat_hydrological_enso %>%
  filter(dry) %>%
  mutate(lat_box = (as.integer(lat) %/% resy) * resy , lng_box = (as.integer(lng) %/% resy) * resy) %>%
  group_by(`Dis No`, lat_box, lng_box) %>%
  summarise(disaster_type = first(`Disaster Type`)) %>%
  group_by(disaster_type, lat_box, lng_box) %>%
  summarise(county = n()) %>%
  filter(county > 0) %>%
  ggplot() +
  geom_sf(data = world_map) +
  geom_tile(mapping = aes(x = lng_box, y = lat_box, fill = county, width = resy)) +
  scale_fill_continuous(type="gradient", low = "white", high = "red") +
  facet_grid(rows = vars(disaster_type))+
  ggtitle("Dry disaster types around the world")

`summarise()` has grouped output by 'Dis No', 'lat_box'. You can override using the `.groups` argument.`summarise()` has grouped output by 'disaster_type', 'lat_box'. You can override using the `.groups` argument.

emdat_hydrological_enso %>%
  filter(dry) %>%
  filter(!is.na(Event)) %>%
  mutate(lat_box = (as.integer(lat) %/% resy) * resy , lng_box = (as.integer(lng) %/% resy) * resy) %>%
  group_by(`Dis No`, Event) %>%
  summarise(disaster_type = first(`Disaster Type`)) %>%
  ggplot() +
  geom_bar(mapping = aes(x = disaster_type, stat = "count")) +
  facet_grid(~Event)

`summarise()` has grouped output by 'Dis No'. You can override using the `.groups` argument.Warning: Ignoring unknown aesthetics: stat

emdat_hydrological_enso %>%
  filter(dry) %>%
  filter(!is.na(Event)) %>%
  mutate(lat_box = (as.integer(lat) %/% resy) * resy , lng_box = (as.integer(lng) %/% resy) * resy) %>%
  group_by(`Dis No`, Event) %>%
  summarise(disaster_type = first(`Disaster Type`)) %>%
  ggplot() +
  geom_bar(mapping = aes(x = Event, stat = "count")) +
  facet_grid(~disaster_type)

`summarise()` has grouped output by 'Dis No'. You can override using the `.groups` argument.Warning: Ignoring unknown aesthetics: stat

Looking at ENSO Value

Create Data as NetCDF

resy <- 2.5

df_to_be_spatialized <- emdat_hydrological_enso %>%
  mutate(lat_box = (as.integer(lat) %/% resy) * resy , lng_box = (as.integer(lng) %/% resy) * resy) %>%
  mutate(month_number = (end_year * 12) + end_month) %>%
  group_by(lat_box, lng_box, month_number, `Dis No`, `Disaster Type`) %>%
  summarize(enso_value = mean(enso_Value), Event = first(Event), num_months_disaster = mean(num_months_disaster)) %>%
  group_by(lat_box, lng_box, month_number, `Disaster Type`) %>%
  summarize(enso_value = mean(enso_value), Event = first(Event), num_disaster_months = sum(num_months_disaster)) %>%
  filter(!is.na(lat_box))

`summarise()` has grouped output by 'lat_box', 'lng_box', 'month_number', 'Dis No'. You can override using the `.groups` argument.`summarise()` has grouped output by 'lat_box', 'lng_box', 'month_number'. You can override using the `.groups` argument.